use 	"./output/cities_expandedpatents_cats", clear

** npi X month-of-year FEs **
egen	idxm = group(gen_id f_m)

** census region X year FEs **
egen	regyr = group(f_yr region)

*=======================*
*= Treatment variables =*
*=======================*

tab 	days_npis if f_myr==tm(1918m9)

gen		daysnpis = days_npis 
gen		lnnpis = ln(days_npis)
gen		longnpis = (days_npis>=90) if !mi(days_npis) // Natural break in data, about 61% have short

foreach v of varlist daysnpis longnpis lnnpis {
	gen 	post_X_`v' 		= `v'*(f_myr>=tm(1918m9))

	gen 	during_X_`v' 	= `v'*( (f_myr>=tm(1918m9)) & (f_myr<=tm(1919m3)) )
	gen 	after_X_`v' 	= `v'*(f_myr>=tm(1919m4))
	gen 	pre_X_`v' 		= `v'*(f_myr<=tm(1917m8))
	
	lab var post_X_`v'		"Post Pandemic $\times$ NPI Length"
	
	lab var pre_X_`v'		"Before Pandemic $\times$ NPI Length"
	lab var during_X_`v' 	"During Pandemic $\times$ NPI Length"
	lab var after_X_`v' 	"After Pandemic $\times$ NPI Length"
}

*=======================*
*= Table 2 Panel A&B: Regression Models =*
*=======================*

** Post ** 

foreach v in daysnpis lnnpis longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_wtd_inv_A `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_A if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_B `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_B if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_C `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_C if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_D `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_D if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_E `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_E if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_F `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_F if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_G `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_G if e(sample)==1
		estadd scalar outmean = `r(mean)'
	eststo: ppmlhdfe pat_wtd_inv_H `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
		sum		pat_wtd_inv_H if e(sample)==1
		estadd scalar outmean = `r(mean)'
	
	estout 	using "$RES/fortable/class_post_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(outmean N, fmt(%9.2f %9.0g) labels("Mean of Dep.\ Variable" "\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	
	capture: ppmlhdfe pat_wtd_inv_Y `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)		
	capture: ppmlhdfe pat_wtd_inv_0 `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)	

	est 	clear
}


foreach v in daysnpis lnnpis longnpis {
	
	local treat pre_X_`v' during_X_`v' after_X_`v'

	eststo: ppmlhdfe pat_wtd_inv_A `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_B `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_C `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_D `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_E `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_F `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_G `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv_H `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/class_duringafterpre_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	
	capture: ppmlhdfe pat_wtd_inv_Y `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)		
	capture: ppmlhdfe pat_wtd_inv_0 `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)	
	
	est 	clear
}
